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Abstract 

Transport properties of three-dimensional self-affine rough fractures are studied by means of an 
effective-medium analysis and numerical simulations using the Lattice-Boltzmann method. The 
numerical results show that the effective-medium approximation predicts the right scaling behavior 
of the permeability and of the velocity fluctuations, in terms of the aperture of the fracture, the 
roughness exponent and the characteristic length of the fracture surfaces, in the limit of small 
separation between surfaces. The permeability of the fractures is also investigated as a function of 
the normal and lateral relative displacements between surfaces, and is shown that it can be bounded 
by the permeability of two-dimensional fractures. The development of channel-like structures in the 
velocity field is also numerically investigated for different relative displacements between surfaces. 
Finally, the dispersion of tracer particles in the velocity field of the fractures is investigated by 
analytic and numerical methods. The asymptotic dominant role of the geometric dispersion, due 
to velocity fluctuations and their spatial correlations, is shown in the limit of very small separation 
between fracture surfaces. 
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I. INTRODUCTION 



The transport of fluids in geological media plays a dominant role in different applications 
such as subsurface hydrology, hydrocarbon recovery and waste storage, and the transport 
properties of such media involves a combination of fluid flow at different length scales. First, 
there is transport at the microscopic level through the pore space of the rock itself. Then, at 
macroscopic length scales, the flow through fractures is, in most cases, the dominant trans- 
port mechanism. Finally, at even larger length scales, the dominant convective transport 
involves the flow through fracture networks. The first case has received considerable atten- 
tion and is relatively well understood and, in fact, it can be placed within the framework of 
transport in porous media [|l], ^ ^. On the other hand, understanding the flow through single 
fractures is clearly a prerequisite to the investigation or modeling of more complex cases such 
as the flow through macroscopic fracture networks. Another key transport process present 
in geological systems, the transport of tracer particles carried by the fluid also requires an 
adequate understanding of the flow properties, in particular, the statistical properties of the 
velocity field. 

Fractures have often been modeled as simple Hele-Shaw cells, with a cubic relation be- 
tween the volumetric flow rate and the average aperture of the fracture, usually focusing 
most of the effort on investigating the flow through a macroscopic network of such fractures. 
However, although a typical fractured rock surface appears fairly smooth, aside from some 
random roughness, suggesting that Poiseuille flow in a straight channel is the appropriate 
model of flow, laboratory experiments |^ and numerical simulations |]^, 0] indicate that 
this classical view of a rock fracture as a straight channel is not adequate to describe the 
fluid flow even at low Reynolds numbers. In fact, more careful analysis showed, in recent 
years, that geological fractures present highly spatially correlated, self-affine surfaces, with 
a roughness exponent ( ~ 0.8 surprisingly constant for different types of rocks whether 



naturally or artificially fractured M, R QG, 11]. In view of these results, theoretical and 



numerical studies introduced the complex geometry of the fractures in order to calculate the 
transport properties of the system. However, the majority of these studies assumed that 
the Reynolds (lubrication) approximation is valid so that the velocity field is given by a 
Poiseuille flow everywhere, with a parabolic velocity profile across the aperture and in the 
direction of the mean flow ||T2| , |13| . 
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In a previous work 0, we showed that, for two-dimensional narrow-fractures, in which 
the roughness amphtude is large compared to the aperture, the Reynolds approximation 
is no longer valid, and proposed a different approach in which the fracture is divided into 
approximately-straight segments with varying orientation angles. This approach allowed 
us to obtain the correction to the flow rate due to surface roughness as a function of the 
aperture and was validated by our numerical simulations. (A similar approach was used by 



Oron and Berkowitz to analyze the case of fractures with contacts between surfaces [Q.) 
In this work, we shall further investigate the case of narrow fractures, extending the results 
to three-dimensional fractures. We shall see that it is possible to obtain analytic expressions 
for the permeability of the fracture and for the velocity fluctuations, in the limit of narrow 
fractures, by means of the effective-medium approximation. We shall also investigate the 
dispersion of tracer particles advected by the flow field within the fracture and its dependence 
on the aperture. 

We shall first, in Sec. II, briefly discuss some basic definitions associated with self-affine 
surfaces and present our model of fractures as the gap between two perfectly matching self- 
affine surfaces. The Lattice-Boltzmann method, used to numerically simulate the fluid flow 
through the fractures is briefly presented in Sec. |T|. Then, in Sec. |IV| we shall investigate 
the permeability by means of the effective-medium approximation and compare the results 
with our numerical results. Finally, in Sec. |V] we investigate the dispersion of tracer particles 
convected by the flow field within the fracture and its dependence on the aperture of the 
fracture. 



II. SELF-AFFINE FRACTURE SURFACES 

Following the work of Mandelbrot [|l^, the application of the fractal model to describe 
surface topography has become widespread as it has been shown in several experiments that 
fracture surfaces, with various materials and fracturing methods, exhibit statistically self- 
affine scaling properties. Particularly relevant to our work are the experiments performed 
with naturally fractured rocks [l^, |19| . This self-affine description of fracture surfaces 



substantially improves the characterization of rough surfaces in terms of surface parameters 
such as the root-mean-square (rms) roughness, rms slope, and density of peaks, in that it 
provides the scaling behavior of these parameters, which are not intrinsic properties of the 
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surface but strongly depend on the sample size 



PI 



A. Scaling properties of self-afRne surfaces 

We shall briefly review here the mathematical characterization of self-affine surfaces, and 
its application to the case of fractured rocks. We consider a rock surface without overhangs, 
as it has been shown that, for natural rock surfaces, overhangs are not prevalent [|T6|, whose 
height is given by a by a single- valued function z{x,y), with the coordinates x and y lying 
in the mean plane of the fracture. 

Self-affine surfaces display scale invariance with different dilation ratios along different 



spatial directions, remaining unchanged under the rescaling 



X Xix y \2y z (1) 

Here we consider disordered media, so these scaling laws apply only in an ensemble or spatial 
average sense. Experiment indicates that for many materials isotropy can be assumed in 
the mean plane, implying that there is only one non-trivial exponent relating the dilation 
ratio in the mean plane to the scaling in the perpendicular direction, i.e., Ai = A2 = A, and 
A3 = A^, where is usually referred to as the roughness or Hurst exponent |jl5[. Therefore, 



the surface height is a homogeneous function, of degree on the mean plane coordinates, 

z{x,y) = \~^z{\x,\y). (2) 

A surface with Hurst exponent close to 1 has high spatial correlation and, as the Hurst 
exponent decreases, heights of nearby points become more uncorrelated, becoming totally 
independent for ^=0.5. In most studies the roughness exponent is found to be close to C = 
0.8, leading to the suggestion that it could be a universal value, independent of the material 
and of the fracture mode ||^. The experimental values obtained with several different types 
of rocks, whether naturally or artificially fractured, are consistent with this suggestion, with 
C lying around C = 0.8 ± 0.05 @, ^ [10|, [Tl|, |l^ (Other values may arise from intergranular 



effects as in the case of sandstone rocks ||2T|.) Therefore, as we shall investigate transport of 
fluids in fractured rocks, we will use a roughness exponent C, = 0.8. Furthermore, it is worth 
noting that the roughness exponent of a granitic fault surface at fleld length scales of several 



meters was found to be = 0.84 [jT8[, thereby supporting the relevance of experimental 
results obtained at laboratory scales. 



B. Characteristic length 



In contrast with the case of self-similar surfaces, where the fractal dimension fully char- 
acterize them, the Hurst exponent is not enough to describe a self-affine surface, and, in 
addition, the amplitude of the fluctuations in the height of the surface is to be specified. 
This amplitude of the fluctuations is usually expressed in terms of the characteristic length 
£, which is the horizontal distance over which fluctuations in height have a rms slope of one. 
Using the self-affine scaling law for the vertical fluctuations in surface height, 

al{r) = {[z{f') - z{f' + f)f) oc r"^ (3) 

and the previously mentioned definition of the characteristic length 

= ^^ (4) 

we can then obtain the variance of the fluctuations in surface height, over any length scale 
r, in terms of the roughness exponent ( and the characteristic length £, 

'r\2C 



crhr 



{r) = {[z{f')~z{f' + f)r)=fQ . (5) 

In principle, the characteristic length i can take any value, however, for fractured rocks 
it is usually found to be very small, and to lie either below the accessible range of length 
scales in the experiment {£ < 0.5/im |jT^) or below the self-affine range of the surface 0, p2| ; 
Therefore, we shall not consider cases in which £ is within the simulated range of length 
scales. 



C. Generation of fracture surfaces 

The self-affine surfaces are generated by a two-dimensional generalization of the Fourier 
synthesis method previously used by us for the generation of self-affine curves 0. Initially, 
an L X L matrix, where each element will represent the height of the discrete version of the 
fracture surface, is generated with statistically independent, Gaussian distributed, random 
numbers. Then, the Fourier transform of this initial random matrix is modulated by a 
power-law high- wave-numbers filter that introduces height-to-height correlations That 
is, if Z{k) is the Fourier transform of the initial random matrix, then the Fourier transform 
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of the matrix representing the fracture surface heights is chosen to be, 



z{k) = k-^-^z{k). 



(6) 



Finally, the whole surface is rescaled in order to get the desired characteristic length i. 
In this work we shall use ( = 0.8 and £ ~ 3 x 10~^. This numerical method generates 
homogeneous and isotropic surfaces and is preferable instead of other numerical methods 



such as the random addition algorithm (for a discussion of this point see Ref. section 
V.A). 




FIG. 1: Example of a numerically generated self-affine surface, with roughness exponent ^ = 0.8, 
characteristic length £ ~ 3 x 10^^, and size L x L = 64 x 64. 

In Fig. 1^ we show one of the numerically generated fracture surfaces. The surfaces 
generated with this method are periodic in x and y and, although the periodicity is not 
physically important it has some computational advantages when calculating the flow field 
using periodic boundary conditions. 

D. Model of fractures 

We shall model fractures as the gap between a numerically generated self-affine surface 
and its replica, which is translated by a fixed distance h in the direction normal to the 
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mean plane of the surface. However, during the fracturing process of the rock, the opposite 
matching surfaces might experience a lateral shift, in addition to the vertical displacement. 
Therefore, we shall also consider this situation by introducing a lateral shift to the upper 
surface either along the direction of the mean flow or perpendicular to it. 

In this work, we will analyze, for both types of fractures, with and without lateral shift, 
the case of narrow fractures, in which the gap size h is small compared to the vertical 
fluctuations of the fracture surface over length scales comparable to the size of the system 
L. In view of Eq. we obtain that fluctuations of the surface height over the whole system 
are given by cr^(L) = iiLjtf and, then the limiting situation of narrow fractures corresponds 
to. 

In our simulations, and for systems size L = 512 the largest gap size that fulfills the 
previous relation is h 

max 

50. 

III. LATTICE-BOLTZMANN METHOD 



The Lattice-Boltzmann Method (LBM) |25| , is particularly suitable to investigate the 
flow of fluids in a highly irregular geometries, and, specifically, to study the various features 
of transport in self-affine rough fractures which are sensitive to their complex geometrical 
structure. In this algorithm, fictitious particles move between neighboring sites on a lattice, 
with suitable collision rules, and the boundary of the flow domain is simply a surface of sites 
where the rule is modified in some way to keep the particles out. We use the Face Centered 
Hypercubic (FCHC-projected) version of the LBM, with a cubic lattice in 3 dimensions and 



19 velocities (D3Q19 in the terminology of |2^). The collision operator is approximated by 



a single relaxation time, the Bhatnagar-Gross-Krook model |^|, and the local equilibrium 
distribution given in is used. This pseudo-equilibrium distribution locally preserves 
mass and momentum values, and is formulated specifically to recover the Navier-Stokes 
equation at large length and time scales. For the no-slip solid boundaries, we use the 
simplest implementation of particle exclusion - the bounce-back rule, in which a particle 
incident on the boundary reverses its direction. Periodic boundary conditions are used for 
the inflow and outflow surfaces. A constant pressure gradient forcing the fluid is added in 
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the X direction, while the gap between surfaces extends over z. 

In what follows all quantities will be rendered dimensionless by the characteristic units 
of the numerical simulation. That is, taking the lattice-spacing as the unit length and the 
simulation time-step as the time unit. Note that, since we are concerned with incompressible 
flows, we do not need to introduce a dimension of mass. The relaxation time is chosen so 
that the kinematic viscosity is = 0.1 in dimensionless units. The pressure gradient is 
Vp = 6.25 10-^, yielding a mean velocity in the range 2.0 10"^ < < 2.0 lO^^ for 
separation between surfaces of the fracture between 8 < /i < 64. Then, for gap sizes h not 
larger than hs — 20, the Reynolds number is Re = < 1 and the flow is governed by 

the Stokes equations, which are invariant under a velocity rcscaling. 

IV. PERMEABILITY OF SELF-AFFINE FRACTURES 

Let us consider first the case in which the two matching surfaces of a fractured rock are 
separated by a distance h in the normal direction, with no lateral shift between the two 
complimentary surfaces of the fracture. In this situation, the aperture of the fracture is 
clearly constant everywhere, but the effective local aperture for fluid flow, i.e., the local 
width of the fracture channel normal to the mean flow direction strongly depends on the 
local angle between the surface and its mean plane. The following theoretical analysis will 
be based on a kind of local lubrication approximation, in which the fracture is divided 
into smaller rectangular blocks. The linear size ^ of the blocks will be estimated as the 
length scale over which fluctuations in the surface height are small compared to the distance 
between surfaces and thus, each of these basic blocks can be considered, approximately, as 
two facing planes separated a distance h and with varying orientation angles with respect 
to the mean plane of the fracture. 

A. Analogue resistor network 

The transport properties of these narrow fractures can be modeled using effective-medium 
ideas. The representation of the fracture in terms of quasi-linear blocks with random ori- 
entation can be mapped onto a regular two-dimensional square lattice, of lattice spacing ^, 
where the disorder only enters the distribution of hydraulic conductances lying on each of 
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the lattice bonds. Moreover, fluid flow is a locally conserved quantity at the lattice nodes 
and therefore, as we shall discuss in this section, this representation of the fracture is com- 
pletely analogous to the classical random resistor networks that model electrical transport 
in disordered media 



Consider the two opposite ends a and 6 of a single quasi-linear block of a fracture, and 
its representation by a conductance joining the corresponding two lattice nodes a and h. 
Moreover, consider that a pressure drop is imposed between the two nodes, and that the 
lattice-bond joining a and h is oriented either parallel or perpendicular to the mean flow. 
Using the Poiseuille result for the flow through this segment of the fracture we have that 

_ [hcos{e)f Apab , . 

12/i [e/cos(^^^ 



where Qa^b is the flow rate from node a to node b, ^ is the linear size of the block, and 6 is the 
orientation angle of the block, with respect to the mean plane of the fracture, in the direction 
of the imposed pressure drop Apab- Let us note, that a slightly different approximation to 
the transport of fluid in a single block can be made by modeling each block as a cylinder 
of diameter h cos{6) and length ^. However, these two approximations will only differ by a 
numerical factor, given that the ratio of the permeability of a cylinder to that of two parallel 
planes is independent of the separation between surfaces, kstraight/kcyiinder = 32/12. 
We can rewrite the previous equation in terms of the bond conductance Qab as 

Apab = -QabQa^b (9) 

where gab is a random variable related to the actual height distribution of the fracture surface 
through the orientation angle 9. 

Furthermore, the fluid is incompressible and flux is conserved at each lattice node, 

Y.Qa^b = Q (10) 

a— >b 

where the sum is over all nodes a connected directly to h. The previous two equations are 
formally equivalent to the equations of an electrical-resistor network with the correspondence 
51 

pressure -v^ voltage (11) 
fluid flux <;=^ electric current 
hydrodynamic conductance -v^ electrical conductance 



Replacing Eq. into Eq. ([T^) we obtain the following system of linear equations for the 
pressure at each node of the network, 

Y,9ab{Pa-Pb)=^ (12) 

Then, the problem is now to solve the previous system of equations in the case where 
the hydrodynamic conductances gab vary according to some probability distribution, and we 
shall discuss now an approximation commonly used to solve this type of problems in the 
physics of disordered media, namely effective-medium theory. 



B. Effective- medium approximation 

The effective-medium approximation (EMA) is a standard method commonly used to 
calculate effective properties of a microscopically disordered medium, in which the random 
microscopic parameters are replaced by a certain mean value, chosen in a self-consistent way 
(a detailed discussion of this approximation can be found eslewhere [^, 

We shall summarize here the main results obtained within EMA for random-resistor 



networks but adapted to flow networks, using the correspondence given in ([llD , and we shall 



make use of these results to compute the permeability of the fracture in Sec. |1VC| and the 
magnitude of the velocity fluctuations in Sec. IV E| . 



—WW 




Sab 
b 



FIG. 2: Representation used to determine the voltage across a conductance gab immersed in a 
uniform network of conductances gm, where gm is the mean EMA conductance. 

Suppose we have an infinite resistor network, with a coordination number z = A and 
individual bond conductances gab, as the one represented in Fig. ^ (the general case where 
z can take any integer value is completely analogous and is discussed elsewhere [|^, W^)- 
The idea then is to choose a mean EMA hydraulic conductance (^m, in a self-consistent way. 
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such that it reproduces the average local field. That is, we require that the extra pressure 
differences induced when one individual conductance reverts from Qm to its original value 
Qab average to zero. In that case, when one bond conductance, oriented along the external 
pressure drop and surrounded by the effective medium, is replaced by an individual random 
conductance having the value Qab, the uniform field solution, in which the pressure drops 
a constant amount Apm across any bond conductance oriented along the imposed pressure 
difference, fails to satisfy fluid flux conservation. In fact, there will be an excess of fluid flux 
in that bond of 

= {gm - gab) Apm (13) 

This excess in fluid flux corresponds to an extra pressure drop Spab, which can be calculated 
in terms of the conductance Gab of the network between points a and b in the absence of gab- 



After some straightforward calculations, the induced pressure drop across ab results |3T 



6pab = Apm ^"" , ^"^ (14) 

gm + gab 

The value of the individual bond conductance gab and hence that of the induced pressure 
difference 6pab is a random variable and, the criterion to choose g^ within EMA is to require 
that the average of 6pab with respect to the probability distribution of local conductances 
gab vanish: 

' gm- g 



(15) 

.gm + g/ 

A detailed study of circumstances under which the EMA is expected to be valid is pre- 
sented in Ref. . 

C. Probability distribution of conductances and mean EMA conductance 

As discussed in the previous section, in order to determine the mean EMA conductance 
gm, we flrst need to obtain the probability distribution of the individual bond conductances 
gab- The conductance of a single block is determined by the angle between the local orienta- 
tion of the surface and the mean plane of the fracture (see Eq. (^). In terms of the height 
difference Z between the two ends of the block, the bond conductance can be written as, 

g{Z) = Goi-^^) , Go = -—- (16) 
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Then, the average over local conductances given in Eq. ([T^) can be performed in terms of 
the distribution of height differences Z, between points of the fracture surface separated 
a distance ^ (note that ^ is the distance between two points lying on the fracture surface 
projected over the mean plane, i.e. the lattice constant.) 

Experimental measurements indicates that the distribution of heights p{Z) can be accu- 
rately described by a Gaussian distribution, at least for low-order moments |T^, Guided 
by these results, the numerical procedure used to generate the fracture surfaces is intended to 



produce self-affine surfaces with Gaussian fluctuations of the height, as discussed in Sec. |1I C 



Moreover, the perturbative analysis we shall present shows that higher moments of the height 
distribution function are related to higher order corrections to Qm and, therefore, the exact 
functional form of the distribution is not relevant but only its low-order moments. 

Finally, assuming that neighboring conductances are uncorrelated, we can rewrite the 
integral equation (^, which defines the EMA conductance gm, in terms of p{Z), 

,„ ,^, g„(i+(z/oy-Go „ 



p(Z) = ^^e.p(-^). (18) 



where 



/27ra2(0 

In terms of the reduced variable x = Z/az{^), which is normally distributed, we obtain. 



where the dimensionless parameter e = (Jz{C)/^ measures the ratio between the average 
magnitude of the fluctuations in the height of the surface cTziC) over a single block, and the 
size of these blocks ^. In general, if the individual blocks of the fracture are smaller than the 
characteristic length ^ ^ £ we have that e < 1, and e > 1 for ^ larger than I. However, since 
we consider the case in which the characteristic length is small compared to the simulated 
length scales {C ~ 10"^ ^ 1, see Sec. |IIB| ), we have that. 



e = 



< 1. (20) 



Therefore, we can then evaluate the mean EMA conductance in a perturbative weak-disorder 
expansion by computing successive terms in the series expansion of gm in the perturbative 
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parameter e, 



12 2 4 

9m — 9m ~l~ 9m^ ~^ 9m^ ~^ 



(21) 



Replacing this series expansion of Qm into Eq. (|T^), and rearranging terms by their order in 
e, we obtain, 



n ) ilm 



9m ~ CIq 
9m + Go 

i9'm + GoV 

2Go 
i9'm + 



[gl + 2{x')gl] 

9m{9m + Go) - {gin) 



+ [2gliGo - g'J] {x') + [gliCo - ^gl)] {x^)] 
+ O(e^) = 0. 



(22) 



Furthermore, since the previous equation must be satisfied for any value of e, then each term 
in the previous equation should be exactly zero, yielding the following expressions for g]^, 



9° 

9m 
9m 



Go, 



-2{x') Go, 



2\2 



[x')+2{x') 



Go- 



(23a) 
(23b) 
(23c) 



Finally, computing the second and fourth moments of the Gaussian distribution given by 
Eq. ( p!8D we obtain the expansion in power series of e for the mean EMA conductance. 



gm = Go [l-2e2 + 5e4 + 0(e6)]. 



(24) 



Let us remark that only the second and fourth moments of the Gaussian probability dis- 
tribution of heights given in Eq. (p!8[) were involved in the previous calculation, and higher 
moments only occur in higher order terms, O(e^). Thus, as already mentioned, this vali- 
dates the use of a Gaussian distribution of heights, which was found to accurately describe 
experimental measurements for low-order moments. On the other hand, to compute higher 
order terms in the e expansion of the electrical conductivity of the network would require 
a different approach, e. g., a perturbative weak-disorder expansion of the conductivity in 
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X 

FIG. 3: Streamlines and front propagation in the gap-averaged flow field. The imposed pressure 
drop and thus the mean flow are along the x direction. The corresponding fracture surface is shown 
in Fig. 1^. The system size is L = 64 and the gap between fracture surfaces is /i = 8. 

terms of the moments of the probability distribution of bond conductances to directly solve 
Eq. (|10|), since EMA would no longer be an accurate approximation ||35| . 



The previous result for the mean EMA conductance gm shows an important feature of 
the fractures, in that the same result, up to second order in e, is obtained if those conduc- 
tances that are perpendicular to the mean flow are eliminated from the resistor network, 
i.e. it predicts a quasi- one- dimensional flow. The two-dimensional character of the network 
only affects higher order terms, O(e^). In fact, deriving from the previous equation the 
permeability of the fracture we obtain. 
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(25) 



which is the same result as obtained in Ref. for two-dimensional fractures (see Eq. (29) 



m 



In Fig. ^ we present the streamlines and the propagation of an initially flat front in the 
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velocity field inside a self-affine fracture, where the flow field was averaged over the gap 
of the fracture. That is, velocities along and perpendicular to the mean flow direction are 
averaged in the z direction over the aperture of the fracture, rendering a two-dimensional 
flow fleld u{x,y). Then, using this flow fleld, the streamlines and the propagation of a flat 
front initially located at x = are computed. The pressure gradient, and therefore the mean 
flow velocity, are along the x direction, {u{x, y)) = Ux x. It can be seen that, as predicted by 
the effective-medium approximation, lateral fluctuations in velocity are small compared to 
the mean flow, and streamlines are approximately straight lines oriented along the direction 
of the imposed pressure difference. However, there is a substantial difference with the 2d 
simulations, in that in the 2d case the gap-averaged mean flow velocity is constant due to flux 
conservation, whereas in the 3d case, as is clear from Fig. |^, fluctuations in the gap-averaged 
velocity occur, and result in the broadening of the initially flat front. 



D. Permeability dependence on the fracture gap 



Thus far we have obtained the relation between the size of the quasi-linear blocks com- 
posing our model fracture and the mean EMA conductance of the system Qm- To further 
obtain the dependence of on the size of the fracture gap, it is flrst necessary to estimate 
the size of the unit blocks ^ in terms of h. The linear size of the unit blocks was deflned as 
the characteristic length over which the channel formed by the two fracture surfaces can be 
considered, in what the flow is concerned, as a straight one. Such an approximation is valid 
when the vertical fluctuations of the fracture surface over the characteristic size of a unit 
block are a small fraction of h, 

^(0 



h 



Q « 1, 



(26) 



where shall be treated as a fltting parameter to our numerical simulations and consitency 
with the previous requirement {C^ <^ 1) will be shown. 
In view of Eq. (Bl) yields. 



h 



Thus, the small parameter e is, in terms of /i. 



\cS 



C-1 



(27) 



(28) 
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FIG. 4: Deviation from the flow rate per unit width in a straight channel qo — q as a function 
of the size of the vertical gap h. We compare results from 2d and 3d simulations and from two 
different sizes L of the system, L = 64 and L = 256. The solid line is a fit to the numerical results 
obtained in the largest 2d system. 

Then, replacing the previous result into Eq. ( pSj ) we obtain the dependence of the perme- 
ability of a three-dimensional self-affine fracture in terms of the separation between surfaces. 



k 



12 



h 



+ 5 



+ 



h 



(29) 



As mentioned earlier, this result is the same, up to second order in e, as the one obtained 
for two-dimensional fractures (see Eq. (33) in |^). In order to test this result numerically, 
we first recast it in terms of the flow rate per unit width q. In a straight channel of height h, 
length L, with fixed pressure drop P, the flow rate per unit width is go = h^P/12fiL. Thus, 
in view of Eq. (p9|), the deviation from go in a fracture of gap size h is given by. 



Qo-q 



]_P (C, 
6/i L 



2-2C 
C 



5C-2 
h — 



(30) 



where is an adjustable parameter expected to be small. 

In Fig. 1^ we compare numerical results, of the deviation from the flow rate per unit 
width in straight channels go, obtained in g — go in 2d and 3d fractures. In both cases, 
the fracture surfaces (fracture curves in the 2d case) have the same roughness exponent 
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and the same characteristic length hence, both have, in average, the same amphtude 
of the vertical fluctuations in surface height. We also show in Fig. ^ the scaling relation 
given by Eq. ( pOf ) with the value of obtained for two-dimensional fractures in Ref. [Q, 
0.1, which as expected is a small number consistent with the approximation of quasi- 
linear blocks discussed before and given by Eq. |2^. It can also be seen that, in agreement 
with the effective medium prediction, both 2d and 3d results are very similar. As discussed 
in Sec. [IV Q the difference between the 2d and 3d cases comes in the fourth order term, 
O(e^), and its correction is of relative magnitude 5e^/2. For the results presented in Fig. ^ 
the largest value of e is e ~ 0.15 (corresponding to the smallest gap size h = 4), hence the 
relative magnitude of the fourth order term is 5e^/2 ~ 0.05, which is consistent with the 
observed similarity between the results obtained in 2d and 3d. 



E. Velocity fluctuations 

In the previous section we showed that the flow rate per unit width in 2d and 3d fractures, 
which have the same characteristic length i, is very similar and, in fact, both have the same 



scaling behavior as a function of the gap size h. However, as was discussed in Sec. [IV C 
there is an important feature of the flow fleld in 3d fractures that is not present in two- 
dimensional simulations, that is the presence of velocity fluctuations even after the fluid 
velocity is averaged over the gap of the fracture. In Fig. ^ we present the fluctuations in 
the gap-averaged velocity, obtained by means of the LBM in a fracture of size L = 64 and 
gap size h = 8. The corresponding fracture surface is shown in Fig. |l] and the streamlines 
corresponding to the gap-averaged two-dimensional velocity fleld are plotted in Fig. ^ 

Let us then investigate the magnitude of these fluctuations in the gap-averaged velocity 
in the direction of the mean flow, i. e. 6ux = Ux{x,y) — U^- Previously, in Sec. [IV B[ , we 
calculated the mean EMA conductance (^m, in a self-consistent way, requiring that the extra 
pressure drop 5pab^ induced between lattice nodes a and h when a bond conductance 
oriented in the direction of the external pressure drop, is replaced by a randomly chosen 
conductance Qah, should average to zero (see Eq. (|15|) ). In a similar way, we can estimate 
the fluctuations in the gap-averaged velocity fleld in the direction of the mean flow, in terms 
of the mean variance of these extra pressure drops. In fact, we might think of the random 
extra-pressure-drops 5pab as the source of the velocity fluctuations. Let us mention, that 
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FIG. 5: Fluctuations in the gap- averaged velocity field The fracture surface used in the numer- 
ical simulation is shown in Fig. and the corresponding streamlines and front propagation are 
presented in Fig. ^. The size of the system is L = 64 and the gap size is h = 8. 

the same result would be obtained if the fluctuations are computed in terms of the mean 
variance of the induced excess of fluid flux from Eq. ([T3|) . 

The variance of the induced extra pressure drops 6pab across a random conductance bond 
Qab oriented in the direction of the flow, is given by, 

1 2> 



9n 



9 



9m + g 



(31) 



Repeating now the procedure we followed to calculate the leading terms of g^, i- e. 
Eqs. ([T7|) , ([18|) and ([19|) , and replacing by its series expansion in powers of e given in 
Eq. (p3|), we obtain 



(32) 



The previous is a general result, in that it is independent of the particular distribution of 
heights of the fracture surface p{x = Z/az{i)). Now, if we replace in this equation the second 
and fourth moments of the distribution by their normal values, (x^) = 1 and (x^) = 3, we 
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FIG. 6: Relative magnitude of the mean velocity fluctuations in the direction of the flow, as a 
function of the gap size h. Two different sizes of the fracture L = 64 and L = 256 are presented. 
Solid line is a fit to the numerical results obtained in the largest system (L = 256). 



obtain, 



(33) 



Finally, using the fact that the relative fluctuations in the pressure drop are equal to the 
relative fluctuation in the velocity, 6p/Apm = 6ux/Ux, we obtain the variance in the gap- 
averaged velocity normalized by its mean value. 



{6ul) {6p^) 



2e' 



(34) 



In order to test this result numerically, let us rewrite the previous equation in terms of 
the gap size h, 



S = CfV2 



h 



(2C-2)/C 



(35) 



where again C/ is an adjustable parameter. 

In Fig. P we show the numerical results obtained for the normalized fluctuations of the 
gap-averaged velocity in the direction of the mean flow, as a function of the distance h 
between unshifted fracture surfaces. We compare results obtained for two different sizes of 
the system L = 64 and L = 256. We find a good agreement with the predicted exponent for 
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FIG. 7: Relative magnitude of the mean velocity fluctuations, both parallel and perpendicular to 
the mean flow, as a function of the vertical separation between surfaces h. The mean flow is in x. 
The size of the system is L = 256. 

the largest system. The fitted exponent is —0.58 ±0.08 and the predicted one is (2 — 2^)/^ = 
—0.5. The adjustable parameter is found to be Cf = 2.1 ± 0.9. Note that can be 
computed from Cf hj means of Eq. obtaining ~ 0.2, which is similar to our previous 
determination and is also consistent with the quasi-linear approximation given by Eq. |2^. 
On the other hand, the scaling behavior of the smaller system deviates from the predicted 
one. We believe that this discrepancy is due to correlations in the velocity fluctuations, and 
we shall show in Sec. |V A| that the correlation length in the velocity fluctuations is, in this 
case, of the order of the size of the system and hence velocity fluctuations are significantly 
reduced. Moreover, we shall also show, in Sec. |V lj|, that as the separation between surfaces 



increases the correlation length also increases, and the mentioned finite-size effect becomes 
more important. This fact explains why the larger deviations from the expected scaling 
behavior observed in Fig. ^ occur at larger separations between surfaces. Let us note, that 
the correlation length of the fluctuations in the gap-averaged velocity is not necessarily equal 
to the size ^ of the quasi-linear blocks. In fact, in the two-dimensional case there are no 
fluctuations in the gap-averaged velocity and still the description in terms of quasi-linear 
blocks was shown to be valid 0. 

In Fig. we compare the fluctuations in the gap-averaged velocity in both directions, 
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along and perpendicular to the mean flow. It can be observed that the scaling of the 
fluctuations in the velocity perpendicular to the mean flow is also given by Eq. with 
a fitted exponent —0.6 ± 0.8. The large uncertainty in the fitted exponent comes from 
the fact that the fluctuations perpendicular to the mean flow have zero mean value and its 
variance presents large variations between different fractures. Finally, it can be observed that 
fluctuations perpendicular to the mean flow are weaker than fluctuations parallel to the flow, 
5ux ~ 5 5uy. The approximately constant ratio between fluctuations along and normal to the 
mean flow, 5ux/5uy ~ 5, is consistent with the conservation of fluid flux. On the other hand, 
the magnitude of the velocity fluctuations are clearly related to the spatial correlations in the 
velocity field, and we will find an analogous asymmetry in the autocorrelation function of the 
velocity fluctuations. Let us mention that similar results are obtained within a macroscopic- 
continuum approach to the problem of transport in heterogeneous porous formations, where 
the magnitude of the fluctuations in the direction of the mean flow are found to be three 



times larger than the perpendicular ones in two-dimensional systems [36 



F. Fractures with shifted surfaces 



Usually, when a rock is fractured its two matching surfaces are laterally shifted, that is 
the displacement between them is not only vertical but also parallel to the mean plane of the 
fracture. We shall now consider this situation, in which the upper surface of the fracture is 
laterally shifted by a vector d = {d\\, dj_) lying in the mean plane of the fracture. However, we 
will only consider the case where the fracture is distinctly open, that is the two surfaces do 
not touch each other at any point. In this the aperture of the fracture is no longer constant 
but becomes a random function of the position ad{x, y) = z{x + d^, y + dy) — z{x, y) + h. Let 
us consider first, the case in which the lateral shift lies in the direction of the mean flow, 
i.e. d = d\\. In Ref. we investigated how such a lateral shift modifies the permeability 
of two-dimensional fractures. We showed that, for large separation between surfaces, such 
that the characteristic size ^ of the quasi-linear blocks is much larger than the lateral shift 
d\\ <^ ^, there is little change in the fracture geometry compared to the unshifted case, 
and that the permeability is asymptotically the same. On the other hand, when surfaces 
nearly touch at some point, the permeability will be dominated by the large pressure drop 
around this point, as the fluid is constrained to flow through this narrow gap. Thus, as the 



21 




h 



FIG. 8: Deviation of the flow rate per unit width from a straight channel, as a function of the 
mean vertical separation between surfaces h. The solid line corresponds to a 2d fracture formed by 
two complimentary surfaces that are simply displaced in the vertical direction, while the dashed 
line corresponds to results obtained when the upper surface of the 2d fracture is also shifted in 
the direction of the mean flow dx = 16. Solid circles correspond to 3d simulations with no lateral 
shift between the matching surfaces. Open triangles and squares correspond to the upper surface 
shifted in the direction of the flow = 16 and perpendicular to it dx_ = 16, respectively. 

surfaces become closer, we found a decrease in the permeability as compared to the unshifted 
case. The case of three-dimensional fractures is somewhat different. For large separation 
between surfaces {d\\ <^ ^; Recall that ^ = ^{h) oc h^^'^) a behavior similar to the unshifted 
case is again expected, since the change in the geometry of the fracture is asymptotically 
negligible, and therefore, the scaling relation given by Eq. (|3D|) should still apply. On the 
hand, when surfaces are close to each other the fluid is no longer forced to flow through the 
narrow gaps, where the minimum separation between surfaces occur, as in the 2d case. In 
three dimension, the fluid can avoid this low permeability regions by flowing around them. 
Thus, we might expect that the fluid rate per unit width q is bounded by the behavior in 
two-dimensional fractures, that is the upper-bound of q given by the two-dimensional flow 
rate in the case without lateral shift, and the lower-bound given by the flow rate per unit 
width in two-dimensional fractures with the same lateral shift, g^l^ < q < qjto- 
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FIG. 9: Schematic representation of the intersection of a small region of the fracture with a plane 
containing both the displacement vector d and the normal to the fracture surfaces. The depicted 
region contains two consecutive quasi-linear blocks of size ^ in the unshifted and shifted cases. The 
local permeabilities are assumed equal in the unshifted case k = k^. 

In Fig. 1^ we present the correction to the flow rate per unit width, go — g, as a function of 
the distance between surfaces, h. In agreement with our previous discussion, the correction 
to the flow rate obtained for three-dimensional fractures lies above the 2d case without 
lateral shift (upper bound for the total flow rate q) and below the two-dimensional results 
obtained for the same lateral shift d = d\\ (lower bound for the total flow rate q). It can 
also be observed that, as we discussed earlier, for large separations between surfaces all the 
results converge to the unshifted situation described by Eq. (pO|). 

In three dimensions, it is also possible to have a displacement perpendicular to the di- 
rection of the flow, d = d±. Let us then investigate how the orientation of the lateral shift 
affects the permeability of the fracture. In Fig. ^we show a schematic representation of the 
intersection of a small region of the fracture, approximated by two consecutive quasi-linear 
blocks of size ^, with the plane that contains both the displacement vector d and the vector 
normal to the mean plane of the fracture. Although extremely simplifled, this schematic 
representation of the fracture shows the effect of the shift on the local permeability. It can 
be seen that, upon a lateral displacement, the unshifted local permeability of a single linear 
block, ko, decreases, or increases, depending on the orientation of the quasi- linear 
block. Speciflcally, an increase (decrease) in the permeability corresponds to < (> 0). 
Furthermore, when the shift is in the direction of the flow the two channels shown in Fig. |^ 
will be in series, i.e. approximately preserving the fluid flux. On the other hand, if the 
shift is in the direction perpendicular to the mean flow the two blocks will be in parallel, 
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FIG. 10: The gap-averaged velocity field 5u{x, y) = u{x, y) — UxX is presented in three different 
cases. On the top is the case with no lateral shift between the surfaces of the fracture. On the 
bottom-left we show the case d = dy = 16, and the case d = d± = 16 is shown on the bottom-right 
corner. In all cases the vertical separation between surfaces is h = 16. 

that is having approximately the same end-to-end pressure drop. It can be shown that, this 
somehow naive model predicts a reduction in the permeability upon a shift in the direction 
of the flow, and a smaller correction in the case where the shift is perpendicular to the flow, 
in agreement with the results presented in Fig. The same effect can be observed in the 
experimental work reported in Ref. [Q, where the flow rate was observed to be larger in 
the direction perpendicular to the shift between the surfaces (see Fig. 4 in [^). It can also 
be observed in the experimental work presented in Ref. that the difference between the 
flow perpendicular and parallel to the shift decreases as the separation between the surfaces 
is increased, which is also in agreement with our results (see Fig. 5 in Let us mention 
that, in our simulations, the largest contrast between q± and q\\, obtained for the smallest 
gap size h = 12, is found to be q±/q\\ ~ 1.24. 



In Fig. 10 we show the effect of the orientation of the shift on the gap-averaged velocity 
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fluctuations, Su{x, y) = u{x, y) — UxX, where the mean flow is subtracted in order to magnify 
the fluctuations in the local velocity. The three different cases presented in Fig. [1^, d = 0, 
d = d\\ and d = d±, have the same vertical gap h, and correspond to the same fracture (same 
self-affine surface). In the case where the lateral shift is perpendicular to the flow d = d±, 
flow channels oriented in the direction of the imposed pressure drop, whereas these oriented 
channels are not present when the shift is along the flow direction d = d\\. This effect is 
sometimes referred to as channeling |Q. 



V. TRACER DISPERSION 



Several tracer-dispersion mechanisms are present in the transport of fluids through porous 
media and the relative importance of these mechanisms strongly depends on the mean 
flow velocity Let us briefly review here the different dispersion mechanisms and their 
dependence on the Peclet number Pe, which measures the ratio between convective and 
diffusive effects and it is deflned, for the fluid flow through fractures, by Pe = hUx/Dm, 
where is the molecular diffusivity, h is the aperture of the fractures and is the mean 
flow velocity. In the case of two-dimensional fractures there are only two contributions to 
tracer dispersion, molecular diffusion, which dominates at very low flow rates, Pe -C 1, and 
is independent of Pe, 0(Pe°), and Taylor dispersion, which is O(Pe^) and therefore becomes 
dominant at high flow rates Pe ^ 1 |]37|, |38|, On the other hand, in three-dimensional 



fractures another mechanism comes into play, that is the presence of spatial fluctuations 
in the velocity fleld. As discussed in the previous section, the effective aperture of the 
fracture is not constant, even in the case when the two complementary surfaces have no 
relative lateral shift. Hence, it is clear that the fluctuations in the local effective aperture 
gives rise to velocity fluctuations, as it was shown in Sec. |1V E| . Moreover, in contrast 
with the two-dimensional case, these velocity fluctuations are present even after the local 
velocity is averaged over the gap of the fracture. One of the main effects of the spatial 
fluctuations in the fluid velocity is that an initially-flat invasion-front of tracer particles will 
become increasingly distorted, resulting in its broadening in time, as it can be observed in 
Fig. ^. This so induced geometric dispersion of tracer particles has been reported in previous 
studies of dispersion in fractures [O, 0| and is completely analogous to that observed in 



three-dimensional porous media |^0[- However, let us note an important difference between 
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previous studies and the present work. In Ref. [|T2|, |T3| the analysis is based on the lubrication 
or Reynolds approximation, where a Poiseuille flow, with a parabolic velocity profile across 
the aperture, is assumed to be locally valid everywhere in the fracture. In this case, there 
are no fluctuations in the local velocity in the absence of fluctuations in the local aperture 
of the fractures, which is in fact the case when there is no lateral shift between surfaces. On 
the other hand, we consider the case of narrow fractures and the lubrication approximation 
is not valid (at least in its simplest version). In fact, in our framework, geometric dispersion 
effects are present even in the absence of any lateral shift between surfaces, due to variations 
in the effective aperture that induce spatial fluctuations in the velocity fleld. This situation 
was investigated in recent experiments, where the broadening and dispersive behavior of 
an invasion-front of tracer particles has been observed even in the case of no lateral shift 
between complementary surfaces 0]. 

In view of our previous discussion, we shall focus on how the fluctuations in the gap- 
averaged velocity affect the dispersion of the tracer particles. The molecular and Taylor 
contributions to the dispersion of tracer particles were discussed in our previous work, in 
the two-dimensional case, and they are not expected to be very sensitive to the fluctuations 
in the gap-averaged flow velocity, in that the molecular diffusion is clearly independent of 
the velocity fleld and the Taylor dispersion is dominated by the gradients in velocity in 
the direction perpendicular to the fracture surface Then, if we only account for the 

geometric contribution to the dispersion of tracer particles the problem can be immensely 
simplifled. In fact, instead of working with the three-dimensional velocity field we can use 
the two-dimensional gap-averaged velocity field u{x,y). The range in which the geometric 
dispersion is the dominant mechanism contributing to the dispersion of tracer particles 
corresponds to intermediate Peclet numbers (intermediate velocities), and the presence of 
such a range of Peclet numbers in self-affine fractures will be discussed in detail at the end 
of this section (see Sec. |VC|) . 
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A. Velocity autocorrelation function 



The dispersion coefficient in the direction of the flow can be expressed in terms of the 
velocity autocorrelation function in time -R||(t) 



/■oo 

Dh = Riiit)dt 
Jo 

= lim / {[u,iO)~U,][u,it)-U,])dt, (36) 



where the integration is over the trajectory of the tracer particles and the average is an 
ensemble average over different realizations of the problem. Furthermore, in Sec. |1V E| we 



showed that the velocity fleld is approximately one dimensional, i.e. lateral fluctuations are 
small compared to the mean velocity {Su/Ux < 0.05). Then, the time autocorrelation func- 
tion of the velocity fluctuations, R\\{t), can be related to the spatial velocity autocorrelation 
function, speciflcally to the marginal spatial autocorrelation function in the direction of the 
mean flow x integrated over the y direction, 

R\\{x) = {[ux{x',y') - Ur,][ur,{x' + x,y') - f4]) , (37) 



by transforming time into position through the mean flow velocity, x = Ux t. In Fig. |TT| we 
show the equivalence between R\\ and R\\, that is R\\{x) = R\\{x/U). Then, we can rewrite 
the dispersion coefficient in terms of the spatial correlation function, 

1 

Dh = - R\\{x)dx. (38) 



U J, 

From the previous equation it is clear that the dispersion coefficient depends on both the 
magnitude of the fluctuations and the length scale of the velocity correlations. Therefore, 
we introduce the characteristic correlation length of the velocity fluctuations, /c, deflned as 

1 /""^ 

lc = — R\\{x)dx, (39) 

Jo 

where au is the rms dispersion in velocity, cr^ = {{ux — UxY)- The velocity correlation 
length Ic measures the typical length over which fluctuations in velocity are correlated and, 
similarly, we can deflne a correlation time Tc = Ic/Ux, which measures the characteristic 
time scale over which fluctuations in velocity remain correlated. Let us remark that Ic is 
not necessarily equal to the previously defined linear size of the quasi-linear blocks ^, a fact 
that becomes clear upon consideration of the 2d case, where one can define a typical size ^ 
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FIG. 11: Spatial velocity-autocorrelation-function R\\{x) and its comparison with R^^{x/Ux), both 
normalized by o"^. The results correspond to a system of size L = 256 with the fracture surfaces 
separated by /i = 8. 



over which the channel formed by the opposing fracture surfaces can be considered straight, 
even though there are no fluctuations in the mean velocity and, therefore, the correlation 
length Ic cannot be defined. 

In an analogous way, we might define the velocity correlation length for the spatial cor- 
relations in the direction perpendicular to the mean flow. 



at 



R_L{y)dy. 



(40) 



However, it can be shown that this integral is exactly zero, due to mass conservation, which 
suggests the presence of negative spatial correlations in the direction perpendicular to the 
flow. 



In Fig. we present both and -R_l for a system of size L = 512 and gap size h = 8. It 
is clear that the spatial correlations in the direction perpendicular to the mean flow decay 
faster than along the flow direction, which is consistent with our previous result concerning 
the magnitude of the fluctuations, in that velocity fluctuations parallel to the mean flow are 
larger than fluctuations perpendicular to it (see Sec. |IVE|) . It can also be observed that both 
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FIG. 12: Spatial velocity-autocorrelation-function, of the velocity component in the direction of 
the flow Ux — Ux, in both the direction of the flow R\\ and perpendicular to it -Rj_. The autocor- 
relation functions are normalized by o"^. The results correspond to a system of size L = 512 and 
fracture surfaces separated by /i = 8. 



correlation functions do not vanish at long distances as it should be in an infinite system. In 
fact, the velocity fluctuations present a positive correlation in the direction of the flow and 
are anti-correlated in the perpendicular direction. This non-vanishing correlation can be 
explained in terms of mass conservation and finite size effects. The gap-averaged velocity Ux 
integrated over a line perpendicular to the mean flow (from y = to ?/ = L) is equal to the 
total flow rate Q divided by the gap size h, and it is a conserved quantity all along the system. 
Therefore, local fluctuations in the velocity Ux should compensate themselves, giving rise to 
negative spatial correlations in Ux along the direction perpendicular to the flow. However, 
this asymptotic negative correlation should decrease as the size of the system increases, as 



it can be observed in our results by comparing Fig. 11, which corresponds to a system of 



size L = 64, with Fig. |T2|, which corresponds to L = 256. This same combination of effects, 
that is mass conservation and the finite size of the system, leads to the positive correlation 
in the direction of the flow observed in Fig. 0. The negative spatial correlation in the 
velocity field along the direction perpendicular to the mean flow, and the fact that the 
velocity fluctuations decay faster in this direction, are also found in the continuum approach 



to transport in heterogeneous porous media [E2|. 
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B. Dependence of dispersion on the gap size 



Let us show first that, as it was assumed in the introduction to Sec. 0, the contribution 
to the dispersion coefficient due to the spatial fluctuations in the velocity is, for vanishingly 
small Reynolds numbers Re, linear in the mean flow velocity Ux- In the Stokes flow approx- 
imation (Re = 0), the flow field becomes independent of the magnitude of the flow rate, Ux 
being just a scaling factor |^3|. Therefore, the dimensionless parameter 6 = cXu/Ux, and the 
correlation length Ic, are independent of Ux- Then, rewriting Eq. (^) we obtain, 

Dh = 5' Ic Ux, (41) 

which explicitly shows the linear dependence of the dispersion coefficient on Ux- Further- 
more, the dispersion of an initially-flat front of tracer particles at a given distance from 
the injection point is independent of Ux- Consider the situation in which a front of tracer 
particles is injected at the inlet section of a fracture (x = 0). The dispersion of the tracer 
front, measured as the mean square displacement of the tracer particles, is then given by, 

{{x - {x)f)^ = 2 Dnt = 26'' IcUxt = (x), , (42) 

and it is clear that the dispersion of the front, at a fixed distance x = (x) from the inlet 
section x = 0, is independent of Ux- In fact, it only depends on = 6' Ic, where Id is 
the dispersion length of the fracture |Q. This fact, that the dispersion of the front is 



independent of the mean flow velocity, was observed in the experiments presented in Ref. 

where it is shown that the front shape depends on the injected volume but not on the 
flow rate, (see Fig. 2 in [H). 

Let us then, investigate the dependence of the dispersion length Id on the gap size h, 
accounting the dependence on h of both the relative magnitude of the fluctuations 6 and 
the correlation length Z^. 

The dependence of 6 on the gape size was discussed in Sec. |1V E| , where we showed that, 



5' oc h-'^'-^^/^. (43) 

On the other hand, in Fig. |1^ we show the correlation l^. as a function of the gap size, 
computed from our numerical simulations using Eq. (|39|), where it can be seen that the 
correlation length increases with the gap size h- Let us mention that, in order to minimize 
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FIG. 13: Correlation length as a function of the gap size h. The size of the system is L = 256. 



finite-size effects, as the previously discussed non-vanishing spatial correlations in the ve- 
locity fluctuations, present in both parallel and transverse directions, the computation of Ic 
was performed in the largest system we could simulate, that is L = 1024. 

The observed decrease in the spatial correlations of the velocity field, as the surfaces 
become closer, might be attributed to a 'screening' mechanism, that is, as h decreases 
the fluctuations in the velocity field become stronger and the velocity tends to decorrelate 
over a shorter distance. An analogous effect occur in porous media flows, where a velocity 
disturbance from a point force decays with a characteristic 'screening' length proportional 
to the square root of the permeability (as seen from the Brinkman equation [^, 0), which 
in our case is proportional to h (the leading order term). 

As discussed earlier, the dependence of on h has two opposite contributions. On one 
hand, the magnitude of the velocity fluctuations, 6, which decreases with h and, on the 
other hand, the correlation length of the velocity fluctuations, Ic, which becomes larger as h 
increases. We found that, as a result of this opposite effects, the dispersion length decreases 
as the gap size is increased, as it can be seen in Fig. 0, where we show the dispersion length 
In as a function of the gap size h. This result is in agreement with the qualitative behavior 
observed in Ref. |Q, where the invasion front of tracer particles becomes smoother as the 
gap of the fracture is increased (see Fig. 5 in |^.) 

Finally, note that underlying the previous discussion is the assumption that the spatial 
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FIG. 14: Dispersion length Id ss a function of the gap size h. The size of the simulated fractures 
is L = 1024. 

correlations in the velocity field decay fast enough so that the integral in Eq. (^) is finite 
and, therefore, the broadening of the tracer front becomes diffusive at length scales larger 
than the correlation length Ic- Analogously, the dispersion of the front is expected to be 
diffusive at time scales larger than the correlation time t^. On the other hand, in previous 
studies of dispersion in self-affine fractures, the slow decay in the spatial correlations of 
the velocity field was proven to induce anomalous dispersion. However, they investigated 
the dispersion of tracer in the lubrication regime, where mean velocity strictly follows the 
fluctuations in the aperture of the fracture and long-range correlations should be expected. 
On the contrary, the lubrication approximation is not valid in the case of narrow fractures, 
and fluctuations in the velocity field are due mainly to the locally random orientation of the 
fracture channel. In Fig. |T3| we present the mean square displacement of an invasion front of 
tracer particles as a function of time. The vertical line shows the correlation time, Tc, after 
which a diffusive behavior should be expected. It can be observed that the initial behavior 
corresponds to the highly correlated motion of the particles and, in fact, the numerical results 
closely follows the solid line that is given by ((Ax)^) = cr^t^. On the other hand, at times 
larger than Tc the velocity begins to decorrelate from its previous values, and the dispersion 
of the front deviates from the quadratic behavior. Moreover, the lower solid line is given by 
a linear fit to the mean square displacement, in the range of times 0.2 < t{Ux/L) < 0.75. 
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FIG. 15: Log-log plot of the mean square displacement of an initially flat front of tracer of particles 
as a function of time. The vertical dashed line corresponds to the dimensionless correlation time 
Tc = Ic/L above which a diffusive behavior is expected. The upper solid line corresponds to 
the initial highly-correlated motion of the tracer particles, (^(Ax)^) = cr^t^- The lower solid line 
corresponds to the best fit of the linear regime, with a dispersion coefficient Dh = (2.7ib0.7) x 10^^. 
(The observed departure from a linear behavior at early times is due to the constant term of the 
fitted linear regime, (A2;)^^q ~ —15. ) The results correspond to simulation in a system with 
L = 1024 and gap size h = 8, and were averaged over 4 different realizations. The time is in units 
of the mean transit time of the medium T = L/Ux- 

The fitted value for the dispersion coefficient is = (2.7 ± 0.7) x 10~^, in agreement 
with the expected value calculated from Eq. (^T]), = (2.0 ± 0.6) x 10^^ However, the 
size of the system is not large enough to observe a large range where the linear, dispersive, 
regime is valid, and therefore the determination of the dispersion coefficient is not accurate. 
Nevertheless, the results allow us to conclude that the mean square displacement grows at a 
much slower rate than in the anomalous regime observed in Refs. |ll2|, [T3[|, where (Ax)^ oc t'^''. 
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C. Tracer dispersion in the three-dimensional velocity field and tracer transit time 
distributions 



The previous analysis of the dispersion problem was based on the assumption that the 
leading contribution to the dispersion of tracer particles comes from the spatial fluctuations 
in the gap-averaged velocity field, which allowed us to map the problem to a two-dimensional 
one. This approximation is only valid for values of the Peclet numbers such that both 
molecular and Taylor dispersion are negligible [0. As we shall discuss, there might be not 
such a range of Peclet numbers, depending on the geometric properties of the fracture. 

The geometric contribution to the dispersion coefficient, given by Eq. (|4l|) , is larger than 
the molecular diffusion term whenever the following inequality holds, 

Dh = 5\U., » ^ Pe » 1 p = 5^^. (44) 

On the other hand, we might expect that Taylor-like dispersion becomes dominant at 
high flow rates, due to the presence of stagnant zones whithin the fracture. In that case, a 
heuristic estimate of the range of Peclet numbers where the geometric contribution generates 
a larger spreading of the tracer front than that induced by the presence of stagnant zones 
is given by, 

Dh = 6\U. » Dt = ^ ^ Pe « /?, (45) 

Then, combining the previous to equations, it is clear that the geometric regime exists only 
for 

P = » 1 (46) 

That is, the product between the magnitude of the fluctuations in velocity and the charac- 
teristic length over which such fluctuations remain correlated should be large. Therefore, 
the existence of such a range of Peclet numbers in which the dispersion of tracer particles 
due to the velocity fluctuations is dominant, would depend on the geometrical properties 
of the fracture. In view of our previous results, in particular the dependence of 6 and Ic 
on h, we might expect that the geometric dispersion would be dominant in the limit of 
narrow fractures, i.e. large fluctuations of the surface height and small separation between 
fracture surfaces. In terms of the small parameter e = CTziO/^^ geometric contribution 
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FIG. 16: Transit time distribution of tracer particles for three different ratios between the transit 
and correlation times, T/tq = 1/2, 1 and5.0, and for three different systems, a Hele-Shaw cell, 2d 
and 3d fractures. The simulation correspond to a system of size L = 512, gap size /i = 16 and a 
transit time distributions are measured at a distance AX = 400 from the injection point. 



will be asymptotically dominant, for any value of the Peclet number, in the limit e ^ 0. On 
the other hand, as the gap of the fracture is increased, the geometric contribution will be 
asymptotically negligible in the limit h ^ oo. 

Finally, let us consider the transit time distribution of tracer particles at high injection 
rates. Previously, we have analyzed the fully-developed dispersion regime, which is valid at 
low injection rates or large fractures. Specifically, let us measure the transit time of tracer 
particles, launched at the inlet section of the fracture, that arrive at the cross-section situated 
at a distance AX. If T is the mean transit time of the tracer particles, T = AX/Ux, and 
T£) is the characteristic correlation time of the tracer velocity, then the Gaussian dispersive 
behavior would be valid for T ^ td- On the other hand, when the injection rate increases 
and the transit time becomes T ^ td, the transit time distribution deviates from a Gaussian 
distributions and exponential tails are generally observed in flow through porous media 
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In the flow through fractures, as well as in the flow in a Hele-Shaw cell, the correlation 
time of the velocity is given by the diffusive time across the gap. For transit times T ^ td 
the tracer particles do not have time to diffuse across the gap and their velocity will remain 
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correlated during their convective motion throughout the system. In this case, stagnant 
or low-velocity zones have the effect of retarding the tracer particles, and give rise to the 
exponential tails, due to the fact that diffusive motion is the only mechanism available for 
the particles to leave these stagnant zones [M] . In Fig. ^ we present the tracer transit time 



distribution for three different ratios between the transit and the correlation times, T/r^j = 
1/2, 1 and 5.0, and for three different systems, Hele-Shaw cells, 2d and 3d fractures. First of 
all, it can be observed that, for T /td ^ 1, all the transit time distributions are Gaussian and 
very similar to each other. On the other hand, as the transit times becomes T/r/j ^ 1, the 
distributions deviate from a Gaussian curve, becoming increasingly asymmetric. It can also 
be observed that both two and three dimensional fractures present slightly more persistent 
tails than in the Hele-Shaw case. The similarity between the distribution in 2d and 3(i 
fractures is in total agreement with our previous results, where the velocity field for three- 
dimensional fractures was shown to be quasi-two-dimensional in the absence of a lateral shift 
between the surfaces of the fracture. On the other hand, we feel that the small difference 
between the transit time distribution in self-affine fractures and in the Hele-Shaw cell might 
be related to the presence of low-velocity zones in the fractures, which are not present in a 
straight channel. This enhancement of the long-time tails due to the presence of low- velocity 
zones should be more evident in the presence of a lateral shift. In Fig. ^ we present the 
transit time distribution of tracer particles, in the same three different cases as in Fig. but 
in this case the upper surface of the fractures is laterally shifted in the direction of the flow 
hy d = dx = 16. First of all, it can be observed that, as in the unshifted case presented in 
Fig. all distributions are Gaussian and very similar to each other for large transit times, 
T /td = 5.0. On the other hand, for much smaller transit times, T/tu = 0.2, a long-time tail 



develops, in particularly in the case of two-dimensional fractures. As shown in Sec. IV F 



two-dimensional fractures present lower permeabilities than the three-dimensional ones in 
the presence of a lateral shift, due to the fact that in the 3d case the fluid can avoid low 
permeability regions by flowing around them. Therefore, the presence of low- velocity zones 
is more important in 2d and thereby the effect of these zones on the long-time behavior of 
the transit time distribution becomes more important. 
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FIG. 17: Transit time distribution of tracer particles for two different ratios between the transit 
and correlation times, T/td = 1/2 and5.0, and for three different systems, a Hele-Shaw cell, 2d and 
3d fractures with opposite surfaces laterally displaced hy d = dx = 16. The simulation correspond 
to a system of size L = 512, gap size h = 16 and a transit time distributions are measured at a 
distance AX = 400 from the injection point. 

VI. SUMMARY AND CONCLUSIONS 

Transport properties of three-dimensional self-afRne rough fractures were investigated by 
means of the effective-medium approximation and numerical simulations using the Lattice- 
Boltzmann method. The numerical simulations verified the scaling behavior predicted by 
the effective-medium approach and, furthermore, allowed us to compute important transport 
parameters of the fractures, such as the dispersion length, and their dependence on the size 
of the aperture. Two different cases were investigated, the unshifted case, in which the two 
matching surfaces of the fracture are displaced in the direction normal to the mean plane, 
and the shifted one, in which the upper surface is laterally displaced either in the direction 
of the flow or perpendicular to it. 

First, we modeled the fractures by a regular two-dimensional square-lattice of bond con- 
ductances, and the lattice spacing and the distribution of bond conductivities were related 
to the geometrical properties of the fracture. Specifically, we related the lattice-spacing to 
the length scale over which fluctuations in the surface height are small compared to the 
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aperture of the fracture, and it was determined in terms of the roughness exponent, the 
characteristic length and the aperture of the fractures. Then, adapting some well-known 
results obtained by means of the effective-medium approximation in the analogous random- 
resistor network, we obtained the permeability of the fracture and its dependence on the 
aperture size in the limit of narrow fractures. We showed that the permeability is, up to 
second order in a perturbative parameter, the same as in two-dimensional fractures. This 
quasi two-dimensional behavior of the transport of fluids through self-affine fractures was 
confirmed by the numerical computations of the streamlines, which presented very small 
lateral fluctuations. A similar behavior was also observed in the experimental work reported 
in Ref. 0], in that the structure developed by the invasion front of tracer particles presents 
very small fluctuations perpendicular to the mean flow in the unshifted case (see Fig. 4a 
in Ref. 0]). Moreover, the scaling behavior of the permeability with the aperture size was 
verifled by our numerical results and, in addition to that, we showed that it is in agreement 
with numerical results performed in two-dimensional fractures. However, we also discussed 
an important difference between the 3d and 2d with cases, namely, the presence of fluctua- 
tions in the gap-averaged fluid velocity in three-dimensional fractures, and we were further 
able to predict the scaling behavior of the velocity fluctuations in the direction of the mean 
flow by means of the effective-medium approximation. The numerical simulations were in 
agreement with this result and, furthermore, showed that the velocity fluctuations in the 
direction perpendicular to the flow have the same scaling but are approximately three times 
smaller in magnitude. Similar results were obtained in the continuum approach to transport 
in heterogeneous porous media l3^, Finally, we investigated the case of shifted surfaces 
and showed that the permeability of the fractures strongly depends on the orientation of 
the shift, which is either in the direction of the imposed pressure drop or perpendicular to 
it, in the limit of narrow fractures. Furthermore, by means of numerical simulations we 
showed that the flow rate per unit width in three-dimensional fractures is bounded by the 
two-dimensional results. Specifically, for a relative shift in the direction perpendicular to 
the mean flow the permeability is slightly affected and lies above the permeability of two- 
dimensional fractures. On the other hand, when the upper surface is shifted in the direction 
of the flow the permeability is largely reduced, but not as much as in the two-dimensional 
case. The latter is due to the fact that in three-dimensional fractures, in contrast with the 
two-dimensional case, the fluid can avoid low-permeability regions by flowing around them. 
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We also presented a simplified representation of a local region of the fracture that, although 
naive in character, captures the effect of the orientation of the shift on the permeability of 
the fracture. The same effect is observed in experiments, in that the flow rate is larger in 
the direction perpendicular to the relative shift between surfaces (see Fig. 4 in [Q). 

In the second part of this work, we investigated the dispersion of tracer particles in self- 
affine fractures. Specifically, we analyzed the dependence of the geometric contribution to 
the dispersion process on the aperture size. First, we simplified the analysis by mapping the 
problem to the dispersion of tracer particles in the two-dimensional gap-averaged velocity 
field. We then distinguised between the two contributions to the dispersion coefficient, 
namely the relative magnitude of the velocity fluctuations and their correlation length. We 



observed that, in agreement with previous studies |36, H2[, the autocorrelation function 



of the velocity fluctuations decays faster in the direction perpendicular to the mean flow. 
Finally, we showed that, even though the correlation length increases with the size of the 
aperture, the dispersion coefficient is asymptotically small in the limit of wide fractures. 
The latter effect is also observed in the experiments presented in Ref. 0] were it was shown 
that the invasion front of tracer particles becomes increasingly smooth as the aperture of 
the fracture is increased. 

Finally, we investigated the dispersion of tracer particles in the fully three-dimensional 
velocity field inside the fractures. Specifically, we discussed the range of Peclet numbers in 
which the geometric contribution to dispersion is expected to be dominant and we showed 
that, depending on the geometric properties of the fracture, there might be no such a range 
of Peclet numbers. However, we found that the geometric dispersion is dominant in the limit 
of narrow fractures, and in general is dominant for large relative fiuctuations of the velocity 
field and correlation lengths larger than the aperture of the fractures. We also investigated 
the transit time distribution of tracer particles, and their dependence on the mean transit 
time. We showed that, as the mean transit time is reduced and it becomes comparable 
or smaller than the correlation time of the tracer velocity, the transit time distribution 
becomes increasingly non-Gaussian, developing long-time tails due to the presence of low- 
velocity zones where the only mechanism for tracer transport is molecular diffusion. In 
general, the transit time distributions are very similar to the case of tracer dispersion in a 
Hele-Shaw cell, except for the case of two-dimensional fractures shifted in the direction of 
the fiow, which present the largest tails probably due to an enhancement of the low-velocity 
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zones by the relative shift between surfaces. In fact, in agreement with the latter results, 
the two-dimensional fractures, with the upper surface shifted in the direction of the flow, 
were shown to have the lowest permeability. 
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